Heavy deformed nuclei in the shell model Monte Carlo method 
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We extend the shell model Monte Carlo approach to heavy deformed nuclei using a new proton- 
neutron formalism. The low excitation energies of such nuclei necessitate calculations at low tem- 
peratures for which a stabilization method is implemented in the canonical ensemble. We apply the 
method to study a well deformed rare-earth nucleus, 162 Dy. The single-particle model space includes 
the 50 — 82 shell plus I/7/2 orbital for protons and the 82 — 126 shell plus O/in/2, 139/2 orbitals for 
neutrons. We show that the spherical shell model reproduces well the rotational character of 162 Dy 
within this model space. We also calculate the level density of 162 Dy and find it to be in excellent 
agreement with the experimental level density, which we extract from several experiments. 

PACS numbers: 21.10.Ma, 21.60.Cs, 21.60.Ka, 27.70+q 



Introduction. The shell model Monte Carlo (SMMC) 
approach [l], Q has been successful in the calculation 
of statistical nuclear properties such as partition func- 
tions and level densities [H, 0|. However, most of the 
SMMC calculations carried out to date have been limited 
to medium-mass nuclei whose deformation is not particu- 
larly large and whose low-temperature properties are well 
described by a single major shell. In such even-even nu- 
clei, the gap to the first excited state is ~ 1 — 2 MeV and 
the ground state can be reached in practice with moder- 
ate values of the inverse temperature (3 ~ 3 MeV -1 . 

For the SMMC to be useful across the table of nuclei, 
it is neessary to show its applicability in heavy nuclei, 
e.g., the rare-earth region, where the deformation in mid- 
shell nuclei can be large and the first excitation energy 
is ~ 100 keV. Such nuclei present a difficult technical 
challenge in SMMC since it is necessary to propagate to 
much large values of j3 ~ 20 MeV -1 . At moderate and 
large values of j3, the propagator becomes ill-conditioned 
and one must stabilize the propagation, keeping its large 
and small scales separated. Stabilization methods were 
developed in strongly correlated electron systems in the 
grand-canonical ensemble [5[. However, nuclear appli- 
cations require use of the canonical ensemble, for which 
stabilization methods are considerably slower. An im- 
portant issue is whether it is possible to describe the 
known rotational behavior of strongly deformed nuclei 
in the framework of a truncated spherical shell model. 
Here we provide an affirmative answer, demonstrating 
our methods for the well-deformed nucleus 162 Dy. This 
is the largest SMMC calculation to date. 

SMMC in proton-neutron formalism. Since protons 
and neutrons occupy different shells, the isospin formal- 
ism is no longer valid and it is necessary to recast SMMC 
in a proton-neutron formalism. A formulation based on 
T z projection was used in Ref. [a|. Here we use a more 
efficient formulation in which protons and neutrons are 
treated explicitly A single-particle orbital a has good 



quantum numbers n, l,j and is 2j + 1-fold degenerate (in 
magnetic quantum number 771) with energy c a . We as- 
sume that the single-particle model space includes Nf 
orbitals for protons (including the magnetic degeneracy) 
and V™ orbitals for neutrons. The two-body interac- 
tion matrix elements are given by Vj P , Vj n and Vj™ for 
proton-proton, neutron-neutron and proton-neutron, re- 
spectively. We first rewrite the two-body interaction in 
a density decomposition by performing a Pandya trans- 
formation for each type of matrix elements to obtain the 
matrices E^, E™ and Ej£ ' . Defining the matrix Ejf as 
the 2x2 block structure with E v ^ and E^ as the diag- 
onal blocks and E p £ , E]P = (Egf as the off-diagonal 
blocks, we have 

a r KM 

where the column vector p KM is composed of both 
proton and neutron densities, and e' a (e' r ) are shifted 
proton (neutron) single-particle energies (the shift orig- 
inates in the Pandya transformation). The matrix 
Ek is real symmetric and can be diagonalized by 
an orthogonal transformation. The quadratic two- 
body term in {T]) can then be written as H' 2 = 
h^Ka^KaY,M(~) M Pk-m{ol)p K m{o), with \ Ka be- 
ing the eigenvalues of Ek ■ The eigenvectors pkm (ck) are 
linear combinations of proton and neutron densities. 

In the Condon-Shortely convention, the time-reversed 
density is given by p K M{ac) = n(-) K+M pK-M iflc), 
where tx = (— ) la+lc is the particle-hole parity. We 
can then rewrite the two-body part of the Hamiltonian 

as H 2 = \ Eifa V KaT,M [QkM^) + R KM( a )] > whcrC 

VfCa — n(—) K \x a , and Qkm, Rkm are the real and 
imaginary parts of pKM (where complex conjugation is 
defined by time reversal). A Hubbard-Stratonovich (HS) 
transformation can be directly applied to this quadratic 
form. The resulting decomposition has a good Monte 
Carlo sign when Vku < for all K, a. The one-body 
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Hamiltonian of the propagator U a in the HS integrand is 
a linear combination of proton and neutron densities and 
the corresponding propagator is a product of a proton 
and a neutron one-body propagators. The computational 
cost is thus smaller than the method used in Ref. Q, in 
which the dimension of the propagator matrix is Nf+N™. 

Stabilization. In SMMC, the evolution operator for a 
given sample is calculated as a product of one-body prop- 
agators of time slice A/3. The number of matrix multi- 
plications increases with (3 and the propagator matrix 
might become ill-conditioned, i.e., the ratio of its largest 
to smallest eigenvalue is too large. Large and small nu- 
merical scales get mixed in the propagation, resulting in 
the loss of important information. 

A method was proposed to stabilize matrix multiplica- 
tion in the grand-canonical formulation [f| . The method 
is based on the decomposition of a matrix M into the 
form M = ADB where A, B are well-behaved under 
multiplication and D is a diagonal matrix whose ele- 
ments are positive numbers containing the various scales. 
In the singular value decomposition (SVD), the matrix 
M has the form M = UDV where U and V are uni- 
tary matrices. In the modified Gram-Schmidt (MGS) 
decomposition M — LDV or M = UDR where L (R) 
is a lower (upper) triangular matrix with diagonal el- 
ements 1. We have adopted the MGS decomposition, 
which can be up to ~ 20 times faster than SVD [8|. In 
SMMC, it is necessary to stabilize the canonical propaga- 
tor. Since the canonical formulation is accomplished by a 
particle-number projection, each term in the quadrature 
sum must be stabilized. 

Choice of model space and interaction. To describe 
the rotational character of a mid-shell rare-earth nu- 
cleus, it is necessary to use a sufficiently large single- 
particle model space. To determine the required single- 
particle orbitals, we consider a Woods-Saxon (WS) plus 
spin-orbit mean-field potential. The spherical orbitals 
of this potential are \ajm) (a represents the remain- 
ing quantum numbers). Introducing an axial deforma- 
tion (3 2 in the WS potential, we determine its eigen- 
states \km) and expand them in the spherical orbitals, 
| fern) = Yl a j c T-aj\ a 3 m )- The spherical occupations 
are then given by r aj = \c% aj \ 2 {n km ), where 

(n>km) are the occupations of the deformed orbitals (1 
below the Fermi energy and above). In a shell model 
approach, we should include in our model space the phys- 
ically important spherical orbitals, while the influence of 
all other orbitals is taken into account by renormalizing 
the interaction. Here we include the orbitals that sat- 
isfy 0.1 < r a j < 0.9 at (3 2 = 0.35. This determined the 
model space to be 0g 7 / 2 , ld 5/2 , ld 3/2 , 2s 1/2 , 0h 11/2 , lf 7 / 2 
for protons, and 0h 11/2 , 0h 9/2 , I/7/2, I/5/2, 2p 3 / 2 , 2p 1 / 2> 
0*13/2: l<?9/2 for neutrons. This model space includes 
orbitals outside the corresponding Ohuj major shells, in 
contrast to Ref. 0- 



As an effective interaction, we use the dominant col- 
lective parts of realistic interactions: monopole pairing 
and multipole-multipole interactions. This interaction is 
similar to the one used in Ref. except that protons and 
neutrons occupy different orbitals 

-53 guP}P»^2xx:{Ox; P +O x . n )-(O x . p +O x . in ): . (2) 

u=p,n X 

Here Pj = E*(-^ +l a^a^ (1/ = p,n) is 
the J = pair creation operator, :: denotes normal or- 
dering, and A; „ = ^ TT E„ i 0'«ll^ a n||i6)[o^ s . w x 

Q>ah;v] 1S a surface-peaked multipole operator [dj m — 
(—iy~~ m aj- m \. We include quadrupole, octupole and 
hexadecupole terms (i.e., A = 2,3,4) with corresponding 
strengths \\ = /£_* k\. The parameter x is determined 
self-consistently [9J and k\ are renormalization factors 
accounting for core polarization effects. 

To determine k 2 , we note that the "slope" of \np(E x ) 
(p(E x ) is the total level density) at higher energies is 
sensitive to \ 2 . We find that a value of k 2 = 2.12 re- 
produces the slope of the experimental \np(E x ) in the 
finite-temperature Hartree-Fock-Bogolyubov (HFB) ap- 
proximation. This value is close to the value of k 2 = 2 
used in Ref. y|. For the octupole and hexadecupole inter- 
actions we take k 3 = 1.5 and k± = 1 

In Ref. 1 we determined the pairing strength to re- 
produce the experimental odd-even mass differences in 
neighboring spherical nuclei using number-projected BCS 
calculations. Following a similar method for spherical nu- 
clei in the mass region Z = 50 — 82, N — 82 — 126, we 
obtain g p = 10.9MeV/Z (g n = 10.9MeV/7V). Here we 
find however that a reduction in the value of g p and g n 
is necessary to reproduce the moment of inertia of the 
ground-state band, and use a reduction factor of 0.77 
(see below). Part of this reduction may be ascribed to 
fluctuations of the pairing fields. 

For the one-body Hamiltonian we use the single- 
particle orbitals of the spherical WS plus spin-orbit po- 
tential. Since the WS potential represents a mean-field 
potential, we determine the bare single-particle energies 
so they reproduce the WS single-particle energies in the 
Hartree-Fock (HF) approximation. 

Ground-state energy and moment of inertia. We 
demonstrate our methods for a typical strongly deformed 
rare-earth nucleus, 162 Dy. To determine the ground-state 
energy it is necessary to extrapolate the thermal energy 
to (3 — 00. We carried out stabilized calculations at 
large (3 values (up to (3 = 20 MeV -1 ) using time slices of 
A/3 = 1/32 MeV- 1 and A^ = 1/64 MeV" 1 . The SMMC 
thermal energy is shown versus T in the top panel of 
Fig. Q] (solid circles). The bottom panel of Fig. [T] shows 
(J 2 ) versus T (J is the total angular momentum). 

A simple way to extract the ground-state energy Eq 
is to assume a ground-state rotational band Ej = Eq + 
h 2 J(J + l)/2X g with a moment of inertia I g . At suf- 
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FIG. 1: Low-temperature thermal energy E (top panel) and 
(J 2 ) (bottom panel) versus temperature T in 162 Dy. The 
SMMC results (solid circles) are fitted to ([3]) (solid lines). 
The dashed lines are results obtained from the lowest five 
experimental bands in 162 Dy. The dashed-dotted line (bottom 
panel) is a fit to the vibrational model result. 



ficiently low temperatures, only the ground band con- 
tributes to thermal observables and a simple calculation 
gives 



E(T) ^E +T 



21 g T/h 2 



(3) 



Fitting a straight line of slope 1 to E(T), we find E n = 
-375.387±0.019 MeV. Comparing with the HFB ground- 
state energy of -Ehfb = —372.263 MeV, we determine a 
correlation energy of E HFB - E = 3.124 ± 0.019 MeV. 

By fitting a straight line 2T g T/h 2 to (J 2 ), we also de- 
termine the moment of inertia T g / ft 2 = 35.8±1.5 MeV -1 
of the ground-state band. This value agrees with the 
experimental value of 37.2 MeV -1 , extracted from the 
excitation energy of the first 2 + state (80.7 keV). 

The SMMC results for (J 2 ) agree with the second re- 
lation in (J3J, derived under the assumption of a rota- 
tional band. This provides evidence that our model space 
is sufficient to reproduce the rotational behavior of this 
strongly deformed nucleus within a truncated shell model 
approach. We also show in Fig. Q] results of a fit to (J 2 ) 
assuming a vibrational model (dotted-dashed line). Our 
SMMC results clearly indicate that the low-lying levels 
of our shell model Hamiltonian are not vibrational. 

To test the validity of the one-band approximation, 
we show in Fig. Q] results for E(T) and (J 2 ) calculated 
using the five lowest experimental bands in 162 Dy (dashed 
lines). The one-band expressions are seen to be valid 
for T < 0.16 MeV down to T » 0.05 MeV. 

Level density. We use the saddle-point expression for 
the level density in terms of the canonical entropy and 
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FIG. 2: Total state density of 162 Dy. The SMMC results 
(solid circles) are compared with the experimental level den- 
sity (dashed line) and the HFB level density (dotted-dashed 
line). Inset: thermal energy E versus /3 (for (3 > 2 MeV -1 ) in 
SMMC (circles) and in HFB (dotted-dashed line). The open 
circles indicate stabilized calculations. 



heat capacity 0], which in turn can be extracted from 
the thermal energy E{j3). Discretization of /3 intro- 
duces systematic errors in E{0) and we found it nec- 
essary to extrapolate to A/3 = using two time slices 
A/3 = 1/32, 1/64 MeV -1 . For /3 < 3.25 MeV -1 we used 
a linear extrapolation in A/3 while for larger values of 
/3 the dependence on A/3 is weaker and we took an av- 
erage value. The results for E(j3) (stabilized for (3 > 3 
MeV -1 ) are shown in the inset of Fig. [5] For comparison 
we also show E(f3) in the HFB approximation (dotted- 
dashed line). 

The SMMC level density is shown by the solid circles 
in Fig. [21 It agrees very well with a composite-formula 
level density (dashed line) , which we extract from several 
experiments (see below). For comparison, we also show 
the HFB level density (dotted-dashed line). We observe 
strong enhancement of the SMMC level density relative 
to the HFB density. Indeed the latter describes only the 
intrinsic states while the SMMC results include all states 
and in particular the rotational bands. 

Experimental Level Density. There are various experi- 
mental data that can be used to determine the level den- 
sity of 162 Dy : an almost complete level scheme at low 

neutron resonance 



and data obtained by the 



excitations (E x < 2 MeV) 1C 
data at E x = 8.196 MeV [13 
so-called Oslo method [13, Qi 

In our level density studies in mid-mass nuclei, we 
used a back-shifted Bethe formula (BBF) to parametrize 
the SMMC level density and compared with similarly 
parametrized data Here we find that at low excita- 
tions a constant temperature formula works better than 
the BBF. We therefore use a composite formula [ItJ that 
combines a constant temperature formula and a BBF 
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FIG. 3: Experimental state density of 1 Dy. Histograms are 
from level counting [To! [Til . [T3 ] , solid squares are renormalized 
Oslo data [l5[ and the triangle is the neutron resonance data 
[ljj. The dashed line is a fit to the composite formula Q. 
Inset (top): the experimental staircase function N(E X ) (his- 
tograms) and its fit to a sixth-order polynomial (solid line). 
Inset (bottom): average level density obtained by a deriva- 
tive of the fit to N(E X ) (solid line) and its fit to a constant 
temperature formula (dashed line). 



energy Em, together with 6 and by, are determined by a 
X 2 fit of \np(E x ) to the Oslo and neutron resonance data. 
We find E M = 1.752 ± 0.036 MeV, which in turn deter- 
mines a = 18.28 ± 0.15 MeV" 1 and A = 0.421 ± 0.014 
MeV. The corresponding composite density is shown by 
the dashed lines in Figs. [3] and [H and is in very good 
agreement with the SMMC state density. 

Conclusion. We have extended the shell model Monte 
Carlo approach to heavy deformed nuclei using a new 
proton-neutron formalism. A stabilization method is im- 
plemented in the canonical ensemble to accurately de- 
scribe the low-energy properties of such nuclei. Applying 
the method to 162 Dy, we show that the spherical shell 
model approach reproduces well the rotational charac- 
ter of this nucleus, as long as a sufficiently large model 
space and an appropriate effective Hamiltonian are used. 
We also calculate the level density of 162 Dy and find it 
to be in excellent agreement with the experimental level 
density, extracted from several experiments. 

This work is supported in part by the U.S. DOE grant 
No. DE-FG-0291-ER-40608 and as Grant-in- Aid for Sci- 
entific Research (B), No. 15340070, by the MEXT, Japan. 
Computational cycles were provided by the NERSC high 
performance computing facility at LBL. 



P(E X ) 



exp [(E x - 

12(B x -A) 5 /4 



El)/Tj] 

,2^/a(E x -A) 



E, r 



< Em 
> Em 



(4) 



The two formulas are matched at an energy Em assuming 
the continuity of the level density and its derivative. 

To determine the state density at low energies we con- 
struct the staircase function N{E X ) (counting number of 
states below E x ) and fit it to a sixth order polynomial 
(top panel of inset to Fig. [3]). Its derivative (solid line 
in botom panel of inset) describes the average state den- 
sity p{E x ). We observe that hxp{E x ) is well fitted by a 
straight line in the range 0.6 < E x < 1.8 MeV, determin- 
ing Ex = -0.73 MeV and Ty = 0.36 MeV. For a given 
Em , the parameters a and A in ([4"]l are determined by the 
matching conditions, s-wave neutron resonance data de- 
termines the sum of the level densities for spin J±l/2 (/ is 
the spin of the target nucleus) at the neutron separation 
energy. Using a spin cutoff parameter of a 2 = XT/H 2 , 
with the rigid-body moment of inertia X 0.015A 5 / 3 ?i 2 

1/2 

and T — [(E x — A) /a] , we obtain the total level den- 
sity at the neutron resonance energy E x = 8.196 MeV, 
where the spin cutoff model is valid [l8( • Additional data 
are from recent gamma-ray spectroscopy experiments [15j 
using the Oslo method [14[, which determines \rxp{E x ) up 
to bo + bxE x with appropriate constants bo and by. The 
measured level density is converted to state density us- 
ing a rigid-body moment of inertia. Since at low excita- 
tions the moment of inertia is reduced from its rigid-body 
value, we only use data for E x > 4 MeV. The matching 
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